require(foreign)

data <- read.dta("panel replication data.dta")
attach(data)

## Estimate model from observed data
summary(lm1 <- lm(chall_increase100 ~ movement + delta_chall_adj_10000 + open1 + extremity + qual_diff + votes1_100 + rep96 + rep98 + dem98 + rep02 + dem02 + rep04 + dem04,data=data))

## generate a vector of simulated coefficients
b.coef = rep(NA, 100000)

## Loop over 100,000 simulations, each time shuffling the values of challenger moderation, estimating the model above, and saving the coefficient for moderation

for(i in 1:100000) {
### randomly "shuffle" the x-values for different y-values
movement.SHUFFLE = sample(data$movement)

### generate new estimates and save coefficients for movement.SHUFFLE
shuffle.lm <- lm(chall_increase100 ~ movement.SHUFFLE + delta_chall_adj_10000 + open1 + extremity + qual_diff + votes1_100 + rep96 + rep98 + dem98 + rep02 + dem02 + rep04 + dem04,data=data)
b.coef[i] <- coef(shuffle.lm)[2]
}

write.dta(as.data.frame(b.coef),"permutation coefficients.dta")

windows(5,5)
plot(density(b.coef),main="Permutation Distribution",xlim=c(-1.5,1.5),xlab="Coefficient Estimates")
segments(-0.36,0,-.36,0.68)





